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ABSTRACT 

We use smoothed particle hydrodynamics simulations of massive protostellar discs to inves- 
tigate the predicted broadening of molecular lines from discs in which self-gravity is the 
dominant source of angular momentum transport. The simulations include radiative transfer, 
and span a range of disc-to-star mass ratios between M^/M* = 0.25 and M^/M* = 1.5. 
Subtracting off the mean azimuthal flow velocity, we compute the distribution of the in-plane 
and perpendicular peculiar velocity due to large scale structure and turbulence induced by 
self-gravity. For the lower mass discs, we show that the characteristic peculiar velocities scale 
with the square root of the effective turbulent viscosity parameter, a, as expected from local 
a-disc theory. The derived velocities are anisotropic, with substantially larger in-plane than 
perpendicular values. As the disc mass is increased, the validity of the a approximation breaks 
down, and this is accompanied by anomalously large in-plane broadening. There is also a high 
variance due to the importance of low-m spiral modes. For low-mass discs, the magnitude of 
in-plane broadening is, to leading order, equal to the predictions from a disc theory and can- 
not constrain the source of turbulence. However, combining our results with prior evaluations 
of turbulent broadening expected in discs where the magnetorotational instability (MRI) is 
active, we argue that self-gravity may be distinguishable from the MRI in these systems if it is 
possible to measure the anisotropy of the peculiar velocity field with disc inclination. Further- 
more, for large mass discs, the dominant contribution of large-scale modes is a distinguishing 
characteristic of self-gravitating turbulence versus MRI driven turbulence. 

Key words: stars: formation, accretion, accretion discs, methods: numerical, radiative trans- 
fer, hydrodynamics 



1 INTRODUCTION 

Protostellar discs are at the heart of both pre-main sequence evo- 
lution and planet formation theory. The accretion of the disc ma- 
terial onto young stellar objects (YSOs) is crucial to explaining 
their phenomenology. This requires angular momentum to be trans- 
ported outwards, so that ma s s may be transp orted inwards (see e.g. 
IPringldl 198 ll lLodatdl2007l ; lArmitagell20 1 ll) . Quantitative models 
of this transport process are needed to describe the evolution of the 
disc's surface density profile, the prot ostar's accretion rate, th e ini- 
tial conditions for planet formation (Chiang & Youdin 201 (J, and 
the subsequent evolution of the planets in the disc. In the simplest 
description, angular momentum transport is assumed to result from 
tur bulence within the disc, l eading to an effective turbulent viscos- 
ity dShakura & Sunvaevll 1973b . 



ac s H. 



(1) 



Here, c 3 is the sound speed, H is the scale height, and a is a param- 
eter that characterizes the efficiency of the trans port. Constraints 
from observed disc lifetimes and accretion rates dHartmann et al.l 
1 19981) are consistent with models in which a ~ 10 

Prostostellar disc angular momentum transport is thought to 
result from two different mechanisms (depending on the loca- 
tion and time in the disc's evolution). In regions of sufficiently 
high levels of ionization, the mag netorotational instability (MRI) 
is a plausible source of turb ulence l lBalbus & Hawlevlll99lL[r998l 
Papal oizou & Neis"or]|2003l ). Further out in the disc (or early on 
when the disc is quite massive), the disc may be cool and massive 
enough tha t the disc's own self-gravity is the dominant source o f 
turbulence dLin & Pringiel ll987ULaughlin & Bodenheimer i 1 19941) . 
This requires that the Toomre Q parameter iToomre 19641 : 
|Purisenetal]|2007l) . 
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where re is the epicyclic frequency (re = Q. for Keplerian discs), and 
E is the disc surface density. Discs with large Q will undergo radia- 
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live cooling towards lower values. The onset of non-axisymmetric 
instability produces spiral structures in the disc, heating it through 
shocks and viscous-scale dissipation. Typically, self-gravitating 
discs will reach a sel f-regulated, qua si-steady state referred to as 
marginal instability (Paczvnski 1978), where the heating and cool- 
ing are in approximate balance, and Q ~ 2. This produces a self- 
sustai ning gravito-turbulence, which can be described by an effec- 
tive a dGammiefcOOlhlLodato & Ricel2 004). alth ough there are cir- 
cumstances where this approximation fails dForgan et al .11201 it) . It 
is also possible that both MRI and self-gravitating turbulence could 
be generated at the same time, particularly in the case of layered ac- 
cretion discs dGammie|l996l : lArmitage et alj200ll ; IZhu et alj2009t 
iMartin & Lubowll20111) . ~ 

Determining observationally which mechanism is driving 
transport at a particular radius in a disc is difficult. Although the 
condition for the onset of self-gravity (Equation [2} is simple and 
well-understood, measuring Q requires an absolute determination 
of the gas surface density as a function of radius. This is hard. 
Most disc mass estimates are derived from dust trace rs, which may 
be systematically biased dAndrews & Williamsi2007r) . Independent 
observational discriminants of turbulence would therefore be very 
valuable. 

One promising route is to detect the turbulent broadening 
of molecular lines in the infrared a nd/or sub-mm. For loca l fluid 
turbulence, elementary arguments dBalbus & Hawlevlll998h sug- 
gest that the characteristic velocity perturbations will be of order 
v ~ \fotC s ~ 0.1c 3 . Using s hearing box simulations of MHD 
turbulence, Si mon et"al] d201 lh show that this is indeed the ap- 
proximate value of velocity perturbations at the disc midplane, but 
that a few scale heights above the midplane they can be as high as 
0.5c s , and that a small fraction of the broadening is in fact super- 
sonic. While extracting turbulent broadening for molecular species 
directly from the simulations would be ideal, doing so requires 
an understanding of the disc's density and temperature fields and 
the abundance of the molecular lines in question. Modelling these 
properties in observations is non-trivial, as is the subsequent mod- 
elling of radiation transport through the disc. Projection and res- 
olution effects must also be considered. Having said this, recent 
observations have successfully constrained turbulent broadening. 
For example, usin g the CO (3-2) trans ition with the Sub Millime- 
tre Array (SMA). [Hughes et al.l d201 lh placed an upper limit of 
v < 0.1 c s in the TW Hydra system and determined a broaden- 
ing value of v ~ 0.4c s in HD 163296. With the advent of more 
sensitive millimetre instruments such as ALMA, these constraints 
will be significantly improved upon. 

In this paper, we present preliminary calculations of turbu- 
lent velocities in global, three-dimensional smoothed particle hy- 
drodynamics simu lations of self-grav itating protostellar discs with 
radiative transfer dForgan et al]|2009h . Our goal is to characterise 
the turbulent velocity distribution to first-order, and identify any 
qualitative features in the turbulent velocity field that might en- 
able the discrimination between alternate mechanisms for driving 
turbulence (although we leave the extraction of synthetic line ob- 
servations for future work). Existing studies of the MRI and of 
self-gravitating turbulence suggest that significant differences may 
exist. Altho ugh large-scale ma gnetic fields play a role in MHD 
turbulence ( Si mon et al.l I2012T) . the MRI can be described to a 
first approximation as a local angular momentum transport mech- 
anism. The sam e is not true of self-gravity in the limit where the 
disc is massive ( Balbu s & Papaloizoull 19991 : lLodato & Ricell2005l : 
Forgar Tet al.l 1201 lh . With this in mind, we seek to identify phe- 
nomenological differences which would allow us to distinguish 



whether line-broadening is MHD-driven or gravity-driven. Such a 
discriminant would be complementary to estimates based upon de- 
rived disc masses, or upo n direct imaging searches for spiral struc- 
ture dCossins et alj|201oh . 



2 METHOD 

2.1 SPH and the Hybrid Radiative Transfer Approximation 

Smoothed Particle Hydr o dynamics (SPH) dLucvl 1 19771 : 
Gingold & Monaghan 1 19771 ; iMonaghari Il992h is a Lagrangian 
formalism that represents a fluid by a distribution of particles. Each 
particle is assigned a mass, position, internal energy and velocity. 
State variables such as density a nd pressure are then ca lculated 
by interpolation (see reviews by Monaghan] 1 19921 1 20051) . In the 
simulations presented here, the gas is modelled using 500,000 SPH 
particles while the star is represented by a point mass particle onto 
which gas particles can acc rete, if they are sufficiently close and 
are bound dBate et al.ll 19951) . 

The SPH cod e used in this w ork is based on the SPH 
code developed by iBate et al I d 19951) which uses individual par- 
ticle timesteps, and individually variable smoothing lengths, hi, 
such that the number of nearest neighbours for each particle is 
50 ± 20. The code uses a hybrid method of approximate ra- 
diative transfer dForgan etai]|2009h . which is built on two pre- 
existing radiativ e algorithms: the poly tropic cooling approxima- 
tion devised bvlStamatellos et alJd2007h. and flux- lim ited diffusion 
(e.g. IWhitehouse & Batel2004l : lMaver et al.l2007l, see lForgan et al.1 
2009 for details). This union allows the effects of both global cool- 
ing and radiative transport to be modelled, without imposing extra 
boundary conditions. 

The opacity and temperature of the gas is calculated using a 
non-trivial equation of state. This accounts for the effects of H2 dis- 
sociation, H° ionisation, He and He + ionisation, ice evaporation, 
dust sublimation, molecular absorp tion, bound-free and free-free 
transitions and electron scatt ering dBell & Lir]| 19941 : iBolev et all 
l2007l : IStamatellos et al.ir2007h . Heating of the disc is achieved by 
P dV work and shock heating. Stellar irradiation is not included as 
a heating source. 

2.2 Initial Disc Conditions 

The gas discs used in this work were initialised with the 500, 000 
SPH particles located between ri n = 10 au and r out = 50 au, 
distributed such that the initial surface density profile is E oc r~ 3 ^ 2 
and the initi al sound speed profi le is c s oc r -1 ' 4 . They were first 
presented in lForgan et al.l d201ll) to evaluate the validity of the a- 
approximation. 

We are primarily interested in considering quasi-steady self- 
gravitating systems, rather than systems that could fragment to 
form bound companions. These initial conditions (in particular the 
small disc radii) were therefore motivated by recent work suggest- 
ing that massive discs will fragment at radii beyond ~ 60 — 70 au 
( Rafik ovll2003: Istamatellos et al.ll2007l: IStamatellos & Whitworthl 
l2008l : IClarkell2009l ; lRice & Armitagell2009l) . This result is consis- 
tent with observa tions that massive disc s tend to have outer radii 
less than 100 au dRodriguez et alj[2005h . A summary of the disc 
parameters investigated can be found in Table Q] Simulation 1, as 
the lowest mass disc in the set, approximates local angular momen- 
tum transport; Simulation 4 exhibits strongly non-local angular mo- 
mentum transport. The strength of non-local effects increases as the 
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Table 1. Summary of the disc parameters investigated in this work. 



Simulation 


Af* (M ) 


M d (M Q ) 


1 


1.0 


0.25 


2 


1.0 


0.5 


3 


1.0 


1.0 


4 


1.0 


1.5 



disc mass is increased, indicated by increasing amplitude in m — 2 
spiral modes (see Figure[T|l. 



2.3 Calculating line of sight velocities 

Observational constraints on turbulence (or, more generally, on de- 
viations from Keplerian rotation) are derived from two-dimensional 
maps of the line of sight velocity field as traced by a particular 
molecular line. Two distinct theoretical quantities that can be de- 
rived from our simulations are relevant for comparing against such 
measurements. One possibility is to compute the density (or vol- 
ume) weighted distribution of line of sight velocity along columns 
that penetrate the disc. This metric would probe only the small- 
scale turbulent velocity field. Our SPH simulations are not best- 
suited to capturing this component accurately, and current observa - 
tional techniques do not resolve these scales l lHughes etai]|201ll) . 
Alternatively, we can subtract off the contribution to the line of 
sight velocity from the mean azimuthal flow, and evaluate the dis- 
tribution of the peculiar velocity field that remains. This peculiar 
velocity field results from the presence of turbulence within the disc 
(since by construction we have removed the velocity of a laminar 
circular disc), but it need not be the same as the local turbulent 
broadening. In particular, for a massive self-gravitating disc, sig- 
nificant contributions to the peculiar velocity field arise from low- 
order azimuthal disc structure, i.e. the spiral arms. We select the 
second metric as it better reflects the observables currently avail- 
able. 

To compute the distribution of line of sight peculiar velocity, 
we raytra ce through the SPH si mulation using the alg orithms de- 
scribed in lAltavetai] J2008h and lForgan & Ricd bold) . Particles 
only contribute to the density and velocity fields along the ray if 
their smoothing volume^ intersect it (see figure[2}. 

As the smoothing volume is spherical, the task of determin- 
ing intersections is reduced to determining the closest approach of 
a ray to a particle (i.e., the impact parameter b, as labelled in fig- 
ure [2}, and comparing it with the particle's smoothing length h. If 
b < 2h, then the particle's smoothing volume is intersected, and its 
properties are used in the calculation. 

Multiple rays are drawn through a region delineated by an an- 
nulus on the disc surface, allowing an azimuthally averaged sam- 
pling of the velocity field while minimising the effects of particle 
disorder and Poisson noise. A total of 250,000 rays are drawn for 
each run, ensuring that almost all particles contained within the 
annulus are intersected by at least one ray. We obtain the mean ve- 
locity in the annulus - in both the vertical direction (v z ), and in the 
plane of the disc (v p ) - and subtract this to obtain the peculiar ve- 
locity field for all particles in the annulus. We then return to the list 



1 The smoothing volume is a sphere of radius 2hi , where hi is the smooth- 
ing length of particle i 




Figure 2. Illustrating the raytracing method. Particles only contribute to the 
estimated line of si ght velocity if their sm oothing volume intersects the ray. 
Figure taken from Forgan & Rice | 2010|) . 

of particles intersected by the 250,000 rays, and perform a density- 
weighted average of the peculiar velocities, in the same fashion as 
Simo n et al. This is a first-order approximation to the emis- 

sivity of the gas (i.e., we expect denser regions to radiate more than 
less dense regions). 

By binning the density-weighted peculiar velocity of each par- 
ticle in the list, we can construct a probability distribution of veloci- 
ties normalised by the local sound speed, P(v/c s ). These linewidth 
probability distributions (LPDs) are not equivalent to any observed 
line profile, but inste ad give the probabil ity of observing a particu- 
lar linewidth. Unlike [Simon et al] ( l201lh . we do not perform time 
averages of these LPDs (although how the LPD varies with time 
can be seen in section [3~3l l. 

2.4 Limitations due to spatial resolution 

The global disc simulations that we use in this work necessarily 
have lower spatial resolution than local, shearing-box simulations. 
Additionally, the use of constant mass SPH resolution elements nat- 
urally concentrates what spatial resolution we do have toward the 
disc midplane, where the bulk of the mass is. These properties of 
the simulations result in two limitations that we need to be mindful 
of. First, we cannot study a radial extent that is too broad, because 
this would result in regions where the low particle density com- 
promises the resolution of physical angular momentum transport 
processes. Second, we cannot study the vertical dependence of the 
turbulent velocity at heights where we have too few resolution ele- 
ments. 

The limitation in the radial extent over which the simulation 
is trust worthy derives from the use of artifical viscosity. While 
required by the SPH code used, the magnitude of artificial viscosity 
must be quantified so that we know where in the disc the artificial 
viscosity is likely to be lower than the effective viscosity generated 
by the gravitational instabili ties. The linear term for the artificial 
viscosity can be expressed a s dArtvmowicz & Lubowll994l ; lMurravl 
ll996tlLodato & Pricel20ld) ; 

fart = 7^ a SPHC s /l, (3) 

where c s is the local sound speed, h is the local SPH smoothing 
length, and qsph is the linear viscosity coefficient used by the 
SPH code (taken to be 0.1). We can define an effective a param- 



© RAS, MNRAS 000,[[H8] 



4 D. Forgan, P. Armitage and J. Simon 



100 



- 1 00 




100 



- 1 00 




Figure 1. Images showing the surface density structure of Simulations 1 (top left), 2 (top right), 3 (bottom left) & 4 (bottom right) after 27 outer rotation 
periods (ORPs). The stellar mass in each case is 1 Mq, and the initial disc masses of 0.25 Mq, 0.5 Mq, 1 Mq and 1.5 Mq respectively. The axis ranges are 
shown in each figure - it is clear that the more massive discs exhibit higher amplitude spiral structures, in particular the m = 2 mode. 



eter associated with the artificial viscos ity by using equation Q3 
dLodato & Ricell2004lForgan et al]2o7lh 



^art 



tCsH, 



and hence combini ng equ ations 

dArtymowicz & Lubovd 19941 : iMurravl 
|201C|) 



(4) 

5J a nd Q gives 



Q art 



Io aspH H- 



(5) 



This shows that where the vertical structure is not well resolved 
(i.e., jj is large), artificial viscosity will dominate. In the simula- 
tions presented he re, this is typically the case inside ~ 10 au (see 
iForgan et al.l201 ll for details), so any data inside this region can not 
be used. We therefore did not initially populate the region inside 10 
au and although particles will move inside 10 au during the course 
of the simulations, we only consider results outside this radius. 

There are als o limitations in the vertical direction. Unlike 
Simo n et al. the simulations are not well resolved above 

two scale heights, so we are unable to comment on how turbulent 
velocities evolve at higher altitudes, and must be satisfied with mid- 
plane data only. Figure [3] shows an attempt to measure the LPD in 



v p at the midplane, and at one scale height above it. The curves 
are remarkably similar, but the poor resolution of higher altitudes 
forbids us from attributing this to any phenomenology of the disc. 



19961 : lLodato & Price] 3 RESULTS 



In this section, we describe the dependence of the discs' peculiar 
velocity fields on the disc mass, radius, and time from the simula- 
tions. 



3.1 Dependence on Disc Mass 

As previously stated, the underlying disc mass may not be correctly 
estimated by observations - can turbulent linewidths help break 
this degeneracy? Figure [4] shows the LPDs for all four simulations 
(where the annulus was placed at 25 au, with rays projected verti- 
cally through the disc). 

From Figure 4 in dForgan et al.l 1201 it) , the time averaged a 
values at 25 au for these four discs range between ~ 0.005 — 0.01. 
Comparing the in-plane (v p ) distributions to these values, we find 
that Simulation 1 has a mode at v jc s ~ y/a, as might be expected 
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Figure 3. The linewidth probability distribution in v p for Simulation 1, with 
r = 25 au. Planar velocities v p are plotted in solid lines, vertical velocities 
v z are plotted in dashed lines. 



Figure 5. Linewidth probability distributions for Simulation 1 , where rays 
are drawn at a variety of radii. Planar disc velocities v p are drawn in solid 
lines, and vertical velocities v z are drawn in dashed lines. 
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Figure 4. Linewidth probability distributions for all four simulations, where 
rays are drawn through r = 25 au. Planar disc velocities v p are drawn in 
solid lines, and vertical velocities v z are drawn in dashed lines. 
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Figure 6. Linewidth probability distributions for Simulation 4, where rays 
are drawn at a variety of radii at zero inclination. Planar disc velocities v p 
are drawn in solid lines, and vertical velocities v z are drawn in dashed lines. 



from a-disc theory. However, increasing the disc mass appears to 
break this relation, which is a symptom of the a-approximation 
itself beginning to break down in the face of non-local angular mo - 
mentum transport (Balbus & Papa loizoulll999l ; lForgan et al]201 ll) . 
The LPDs change shape dramatically, developing strong peaks and 
strongly asymmetric profiles. These profiles appear to be tracing 
the low-order m — 2 spiral arms as can be seen in the lower panels 
ofFigureQ] Despite this, there is a strong cutoff at v/c s — 1. This is 
due to the subtraction of the average planar motion, which for self- 
gravi tating discs at marginal instability is v /c s = 1 dCossins et al.l 
120091) . 

The v z distributions are much more orderly, with the mode 
increasing by around one order of magnitude as the disc mass is 
increased by a factor of 5. The shape of these profiles is also rea- 
sonably constant. 

3.2 Dependence on Disc Radius 

Having seen the striking change in LPD as the disc mass is in- 
creased and low-m spiral modes begin to significantly affect line 
broadening, it is instructive to compare how LPDs change as dif- 
ferent disc radii are probed. 



Figure [5] shows LPDs for Simulation 1, drawn at a variety of disc 
radiiQ In this low mass case, the difference in peak between planar 
and vertical LPDs remains constant across all radii. Indeed, both 
LPDs remain very similar across all radii. The a profile of Simula- 
tion 1 begins to flatt en out towards a constant value beyond r ~ 30 
au (see iForgan et alJlibT l|), which would presumably explain why 
the peaks of the LPDs tend towards similar values (as the charac- 
teristic velocity v/c s ~ \/ot). 

The situation is somewhat different for Simulation 4 (Figure[6]l. As 
the radius increases, the presence or absence of a spiral arm inside 
the annulus strongly affects the observed LPD for v p . Annuli con- 
taining an arm have a strong spike, with a much narrower dispersion 
in v/c s . Annuli which do not contain arms possess broader LPDs. 
Notably, the v z distributions remain unchanged, although this may 
be symptomatic of low vertical resolution, and as such we should 
be wary of denoting this as a phenomenological effect. 



2 Due to numerical limitations, w e are restricted to a relatively small range 
in radii compared to the studies of l Simon et all 1201 ll) ; see Section l2~4l 
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Figure 7. Linewidth probability distributions as a function of time for Sim- 
ulation 1. Rays are drawn at r = 25 au. Planar disc velocities v p are drawn 
in solid lines, and vertical velocities v z are drawn in dashed lines. 

3.3 Time Dependence 

Finally, it is expected that linewidths should be able to probe the 
inherent variability in these self-gravitating discs. As is shown in 
iForgan et al.l d20 1 IT) , Simulation 1 exhibits low variability, and Sim- 
ulation 4 exhibits high variability. These features can be detected by 
calculating the LPD at a series of timesteps throughout the disc's 
evolution. 

The LPDs of Simulation 1 remain stable over several thousand 
years (Figure|7](, with the differences for both v p and v z being lit- 
tle more than extremely small shifts in the value of the peak. The 
systematic separation between v p and v z is maintained as a result. 

By contrast, the v p LPD for Simulation 4 (Figure[8j fluctuates 
much more noticeably, with both the shape of the LPD and the 
amplitude of the peaks shifting over time. This is again due to the 
passage of strong m = 2 features through this particular radius. 

Observers hoping to construct LPDs by observing objects 
should take note. These results would suggest that such efforts 
would be successful for low mass discs, but the increased variabil- 
ity would rule this method out for high mass discs. Conversely, the 
simulations show that variability can be seen in higher mass discs, 
even at this comparatively low resolution. 

Again, we see very little change in the v z LPD across all disc 
masses, but we must be cautious as these distributions are the most 
sensitive to resolution issues. 

3.4 Dependence on Numerical Resolution 

Finally, we compare the effects of increasing particle number on 
the resulting LPDs. Simulation 1 was rerun with 1 million parti- 
cles (double the previous value), and LPDs were investigated at the 
same instant as a function of radius. Due to computing constraints, 
the high resolution simulation was not run for the same number of 
outer rotation periods as the low resolution simulation. In Figure 
[9] we compare the two resolutions at the same simulation time of 
approximately 3500 years. 

The general features of the LPDs are robust - the position of 
the peaks remains at similar values. The v z distributions remain 
very much the same, despite increasing the vertical resolution by 
around 25%. Strong differences can be seen in the outer disc radii 
(blue and green curves). This is due to extra resolution having a 
stronger effect at lower particle densities. The inner disc regions 



0.1 2 




10.0000 



Figure 8. Linewidth probability distributions as a function of time for Sim- 
ulation 4. Rays are drawn at r = 25 au. Planar disc velocities v p are drawn 
in solid lines, and vertical velocities v z are drawn in dashed lines. 

still display systematic separation of v p and v z , although the am- 
plitude of the peaks in v p do change by around 10%. 



4 DISCUSSION & CONCLUSIONS 

We have presented preliminary studies of turbulent line broadening 
due to self-gravity using global smoothed particle hydrodynamics 
(SPH) simulations of protostellar discs. With radiative transfer in- 
cluded in these calculations, we carried out raytracing techniques to 
measure the local velocity components in the disc plane and trans- 
verse to the disc plane (relative to the local sound speed). From 
these raytracing measurements, we construct a linewidth probabil- 
ity distribution (LPD) which gives the probability of detecting lines 
broadened to a given width. These simulations cannot resolve tur- 
bulence at scales below the smoothing length, and artificial viscos- 
ity effects will also smooth out turbulence to a lesser extent. That 
being said, we have detected some promising diagnostics that can 
be compared to current MRI-driven turbulence simulations but that 
will require follow-up with high resolution studies. 

First, we have found that the LPD is noticeably mass depen- 
dent, with the maximum value of v/c 3 increasing towards v/c s — 
1, where there is a strong cutoff. Non-local angular momentum 
transport is evident from strong peaks in the distribution in v p being 
caused by low-m spiral structure impinging on the disc radii being 
observed. This behaviour becomes apparent when the same disc is 
compared at various radii of observation. 

Low mass discs which possess weaker, high m spiral struc- 
ture retain similar LPDs at all radii, with a systematic separa- 
tion between the peak of the v z distribution and the peak of the 
v p distribution; v p peaks near ~ 0.04e s , whereas v z peaks near 
~ 0.005 - 0.01c s . 

The results from this work can be compared with the work 
of ISimon et alj d201 ll) as a way to extract potential differences be- 
tween the LPDs of MRI-driven turbulence and motions due to self- 
gravity. First, comparing the actual mid-plane values of v/c s , we 
find roughly the same (though somewhat lower) peak v p /c s in our 
low m ass discs compar ed to MRI-driven discs in the ideal MHD 
limit dSimon et al.l201 lb . However, the peak values of v p /c a in this 
work are slightly larger than the dead-zone values of the 4 au sim- 
ulation of Simon et al. 

The largest difference between the LPDs of low mass self- 
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Figure 9. Linewidth probability distributions for Simulation 1 at different numerical resolutions, where rays are drawn at a variety of radii. Planar disc 
velocities v p are drawn in solid lines, and vertical velocities v z are drawn in dashed lines. The left panel indicates the 500,000 particle run, the right panel the 
one million particle run. 



gravitating discs and MRI-driven turbulent discs is the significant 
separation of v z and v p in the self-gravitating case; this robust re- 
sult is not present in the MRI-driven simula tions as v z and v p have 
roughly the same LPD dSimon et al.ll201 it) . As such, it would ap- 
pear to be a good diagnostic for distinguishing between these two 
angular momentum transport mechanisms in the limit of low disc 
mass. 

In the high disc mass limit, the LPD fluctuates strongly in v p 
as the observation radius is changed because spiral arms move in 
and out of the annulus. No such fluctuation is detected in v z , but 
this is possibly due to resolution effects smoothing them out. The 
high time variability exhibited by these arms is also apparent in 
the LPDs as measured across time. This is another feat ure that was 
not observed in the simulations of Isimon et al.l d201ll) . The angu- 
lar momentum transport in these simulations was highly local and 
there were no low m features observed in the LPDs. However, these 
particular simulations were carried out in the local limit; going to 
larger scales by performing global MRI simu lations may introduce 
the presence of low m features into the LPD (Beckwit h et alj201 ll ; 
ISimon et al]|2012h . 

In conclusion, the global simulations carried out in this work 
suggest that one could potentially distinguish between MRI-driven 
turbulence and self-gravitational effects by observing the varia- 
tion of turbulent line width with radius (to potentially probe low 
m structures) and inclination angle. We must reiterate, however, 
that our simulations are quite low in resolution, and as such, im- 
portant phenomenology may not be exhibited. In addition to this, 
stellar irradiation - an important driver of line emission from the 
disc surface - is not included as a radiation source. We recommend 
that local simulatio ns be carried out in much the same manner as 
Simo n et al.l d201 lh . These simulations would be capable of con- 
firming the features of gravito-turbulence seen in the low-mass 
discs, whose angul ar momentum transport is locally determined 
(Forgan et al. 2011). This includes the systematic separation of v z 
and v p . Also, the increase in resolution afforded by local simula- 
tions would allow the study of gravito-turbulent broadening as a 
function of disc altitude, giving a further means of comparison to 
MRI turbulence. Stellar irradiation must also be considered in fu- 
ture simulations used for linewidth studies.Finally, it is clear that 
more sophisticated attempts to make synthetic observations of the 
linewidth broadening are required to determine the detectability of 
these features using current and future telescopes. 



It is also worth noting that while theory indicates that self- 
gravitating protostellar discs should exist only in the early embed- 
ded phase of protostar formation, this does not present an insur- 
mountable obstacle to observing these features. ludicious selection 
of a line with high critical density for activation should allow the 
midplane to be probed without contamination from the envelope. 
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